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We calculate the secular changes of the orbital parameters of a point particle orbiting a 
Kerr black hole, due to the gravitational radiation reaction. For this purpose, we use the 
post-Newtonian (PN) approximation in the first order black hole perturbation theory, 
with the expansion with respect to the orbital eccentricity. In this work, the calculation 
is done up to the fourth post-Newtonian (4PN) order and to the sixth order of the 
eccentricity, including the effect of the absorption of gravitational waves by the black 
hole. We confirm that, in the Kerr case, the effect of the absorption appears at the 
2.5PN order beyond the leading order in the secular change of the particle’s energy and 
may induce a superradiance, as known previously for circular orbits. In addition, we find 
that the superradiance may be suppressed when the orbital plane inclines with respect 
to the equatorial plane of the central black hole. We also investigate the accuracy of the 
4PN formulae by comparing to numerical results. If we require that the relative errors 
in the 4PN formulae are less than 10“®, the parameter region to satisfy the condition 
will be p > 50 for e = 0.1, p > 80 for e = 0.4, and p > 120 for e = 0.7 almost irrespective 
of the inclination angle nor the spin of the black hole, where p and e are the semi-latus 
rectum and the eccentricity of the orbit. The region can further be extended using an 
exponential resummation method to p > 40 for e = 0.1, p > 60 for e = 0.4, and p > 100 
for e = 0.7. Although we still need the higher order calculations of the PN approximation 
and the expansion with respect to the orbital eccentricity to apply for data analysis of 
gravitational waves, the results in this paper would be an important improvement from 
the previous work at the 2.5PN order, especially for large p region. 

Subject Index EOl, E02, E31, E36 


1. Introduction 

The gravitational two-body problem is a fundamental issue in general relativity. This also 
attracts great interest in gravitational wave physics because binary inspirals are promising 
sources of gravitational waves which are expected to be detected directly by ongoing grav¬ 
itational wave observatories in the world. Understanding the dynamics of binary system is 
required to predict the emitted gravitational waveforms accurately for efficient searches of 
the signal in observed data. 

One of major approaches for this purpose is the gravitational self-force (GSF) picture in 
the black hole perturbation theory. In this picture, a binary is regarded as a point mass 
orbiting a black hole and the dynamics can be described by the equation of motion of the 
mass including the effect of the interaction with the self-field, that is, the GSF. After the 
formal expression of the GSF was presented by Mino, Sasaki and Tanaka [1] and Quinn and 
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Wald [2], a lot of efforts have been devoted to develop practical formulations and methods to 
calculate the GSF (for example, refer to [3] for the formulation of GSF, [4, 5] for the recent 
progress in practical calculations of GSF). 

Although a lot of progress has been made, however, it is still challenging to calculate 
the GSF directly for general orbits, especially in Kerr spacetime. Practical calculations of 
the GSF with high accuracy will require a huge amount of time and computer resources 
mainly because of the regularization problem induced by the point mass limit. Therefore it 
is important to develop a way to reduce the cost of computing the GSF. The two-timescale 
expansion method [6] gives a hint for it: assuming that a point mass does not encounter 
any transient resonances {e.g. shown in [7]), the orbital phase, which is the most important 
information to predict the waveform, can be expressed in the expansion with respect to the 
mass ratio, r], as 




$( 0 ) + 0 (^ 2 ) ^ 


( 1 ) 


where and are quantities of order unity. The leading term, can be calculated 
from the knowledge up to the time-averaged dissipative piece of the first order GSF, corre¬ 
sponding to the secular growth. The calculation of this secular contribution can be simplified 
significantly by using the radiative field defined as half the retarded solution minus half the 
advanced solution for the equation of the gravitational perturbation [8-10], i.e. the adia¬ 
batic approximation method, because the radiative field is the homogeneous solution free 
from the divergence induced by the point mass limit. This method allows us to calculate 
the leading term accurately without spending huge computational resources. On the other 
hand, the calculation of requires the rest of the first order GSF (the oscillatory part 
of the dissipative GSF and the conservative GSF) and the time-averaged dissipative piece 
of the second order GSF. There is no simplification in calculating these post-1 adiabatic 
pieces at present. Since is subleading, however, the requirement of the accuracy is not 
so high compared to that of the leading term. This fact suggests that it is possible to reduce 
the computational cost by using a suitable method with an appropriate error tolerance to 
calculate each piece of the GSF (for example, a hybrid approach is proposed in [11]). 

In this work, we focus on the time-averaged dissipative part of the first order GSF, which 
has the dominant contribution to the evolution of inspirals, and present the analytic post- 
Newtonian (PN) formulae. So far, several works in this direction had been done for two 
restricted classes of orbits: circular orbits and equatorial orbits. (See [12] and references 
therein for early works in 1990’s). Recently, thanks to the progress of computer technology, 
the very higher order post-Newtonian calculations can be possible for circular equatorial 
orbits: the 22PN calculation of the energy flux is demonstrated in Schwarzschild case [13], 
and the IIPN calculation in Kerr case [14]. There is also the calculation of the secular GSF 
effects for slightly eccentric and slightly inclined (non-equatorial) orbits [10], and later it 
had been extended to orbits with arbitrary inclination [15], where the PN formulae of the 
secular GSF effects are presented in the expansion with respect to the orbital eccentricity. 
However, the calculation in [15] has been done only up to the 2.5PN order with the second 
order correction of the eccentricity. Also the absorption to the black hole is ignored there. 
The main purpose of this work is to update the results in [15] up to the 4PN order and the 
sixth order correction of the eccentricity, including the effect of the absorption to the black 
hole. 
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This paper is organized as follows. In Sec.2, we give a brief review of the geodesic motion of 
a point particle in Kerr spacetime, the gravitational pertnrbations induced by the particle, 
and the adiabatic approximation method of calcnlating the secular effect of the GSF. In 
Sec.3.1, we present the PN formulae of the secular changes of the the energy, azimuthal 
angular momentum, and Carter parameter of the particle due to the gravitational radiation 
reaction in the expansion with respect to the orbital eccentricity. In Sec.3.2, we investigate 
the accuracy of our PN formulae by comparing to numerical results given by the method 
in [16-18], which can give each modal flux at the accuracy about 14 significant figures. In 
Sec.3.3, we implement a resummation method to the PN formulae given in this work in order 
to improve the accuracy. In Sec.3.4, we discuss the convergence of the analytic formulae as 
the PN expansion and the expansion with respect to the eccentricity. Finally, we summarize 
the paper in Sec.4. For the readability of the main text, we present the PN formulae for the 
orbital parameters, the fundamental frequencies, the orbital motion in Appendices A and B, 
which are used in calculating the secular changes of the orbital parameters. We also present 
the PN formulae for the secular changes of an alternative set of the orbital parameters in 
Appendix C. Throughout this paper we use metric signature (—h ++) and “geometrized” 
units with c = G = 1. 


2. Review of formulation: Adiabatic radiation reaction 

The orbital evolution of a point particle due to the time-averaged dissipative part of the 
GSF is often described in terms of the secular changes of the orbital parameters. In order to 
calculate the changes, we need the information on the first order gravitational perturbations 
induced by the particle when it moves along the background geodesics. In this section, 
we review the geodesic dynamics of a point particle in Kerr spacetime, the gravitational 
perturbations induced by the particle, and the adiabatic evolution of the orbital parameters. 


2.1. Geodesic motion 

The Kerr metric in the Boyer-Lindquist coordinates, {t,r,9,(p), is given by 


/ 2Mr\ 

g^.dx^dx'^ = 


dt^ — 


AMar sin^ 9 


-dtdif + -^dm 


+{r^ + a^ + 


2Ma?r 


sin^ 9 ) sin^ 9dLp ^, 


( 2 ) 

where S = cos^ 0, A = r^ —2Mr + a^, M and aM are the mass and angular 

momentum of the black hole, respectively. 

There are two Killing vectors related to the stationarity and axisymmetry of Kerr space- 
time, which are expressed as = (1,0,0,0) and = (0,0,0,1). In addition, it is known 
that Kerr spacetime possesses a Killing tensor, where P and are 

the Kinnersley’s null vectors given by 


:= 


-I- a? 


A ’^’°’A 


:= 


A a 
2E 


(3) 


For the geodesic motion of a particle in Kerr geometry, there are three constants of motion 
related to the symmetries: 


E := L := Q := 


(4) 
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where u" is the four velocity of the particle. E and L correspond to the specihc energy and 
azimuthal angular momentum of the particle respectively. Q is called as the Carter constant, 
which corresponds to the square of the specihc total angular momentum in Schwarzschild 
case. These specihc variables are measured in units of the particle’s mass, One can recover 
the expressions in the standard units as 

E := ixE, L := /iL, Q := fj?Q. (5) 

There is another dehnition of the Carter constant, C = Q — (aE — L)^, which vanishes for 
equatorial orbits. In this paper, we use C as one of the orbital parameters, instead of Q. 

By using these constants of motion, the geodesic equations can be expressed in the following 
form as 

f dr , (dcos6\^ 

U) (—) =e(cos«), (6) 

^ = V,r(r)+ ¥„(«), ^ = VM + VW(«). (7) 

where we introduced a new parameter A through the relation dX = dr/S, and some functions 
as 


P{r) 

:= E{r^ + a^) — at, 

(8) 

R{r) 

:= [P{r)f - + {aE - Lf + C], 

(9) 

0(cos 9) 

:= C-{C + 0^(1 - E^) + L^) cos^ 9 + a^{l - E^) cos^ 9, 

(10) 

Vtrir) 

2 2 

:= ^ P{r), Vte{9) := a{aEs\T?9 L), 

(11) 

Vprir) 

:= IPir), V 

(12) 


A generic geodesic orbit in Kerr spacetime can be characterized by three parameters, 
{E,L,C} In the case of a bound orbit, we can use an alternative set of parameters, 
{rp,ra,6mm}, instead of {E,L,C}, where and are the values of r at the periapsis and 
apoapsis and dmin is the minimal value of 6, respectively. Using this set of parameters, we can 
describe the range in which the motion takes place as < r < and dmin < d < tt — dmin- 
There is another useful choice of parameters used in [19], {p, e, i}, defined by 

2rpra Ta-Vp L 

V '■= WT7-e := -cosi := , =. 

M{ra + rp) Ta + Tp y/L^ + C 

By analogy to the parametrization used in celestial mechanics, p, e, i are referred as 
semi-latus rectum, orbital eccentricity, orbital inclination angle, respectively. For later con¬ 
venience, we also introduce Y = cos i and v = \/ljp- Since v corresponds to the magnitude 
of the orbital velocity, it can be used as the post-Newtonian parameter. For example, we 



^ Strictly speaking, the orbit is also characterized by the initial position of the particle. However, 
the time-averaged dissipative part of the first order GSF does not affect the initial position (the 
other parts of the first order GSF and the higher order GSF will do) [6]. Also the secular changes 
of {E,L,C} does not depend on the initial position. Hence we do not need the information on the 
initial position to describe the secular evolution of the orbit at the order considered in this paper. 
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call the 0(?;®)-correction from the leading order as the fourth order post-Newtonian (4PN) 
correction. 

It is worth noting that, by introducing A, the radial and longitudinal equations of motion in 
Eq.(6), are completely decoupled. For an bound orbit, therefore, the radial and longitudinal 
motions are periodic with the periods, {A^, A^}, defined by 


hr = 2 



dr 

TWy 


* . do 

Ae = 4 / — 

vW) 


This means that these motions can be expressed in terms of Fourier series as 

OO 

r{X) = pYI COS 

nr=0 

OO 

cos0(A) = \/\ — (5ne sinneOgA, 

ng=0 


(14) 


(15) 

(16) 


where and fig are the radial and longitudinal frequencies given by 




27r 




2-k 

A^’ 


(17) 


and we choose the initial values so that r(A = 0) = and 6{X = 0) = 7r/2. ^ 

Since the temporal and azimuthal equations of motion in Fq.(7) are divided into the r- 
and 0-dependent parts, the solutions can be divided into three parts: the linear term with 
respect to A, the oscillatory part with period of A^, and the oscillatory part with period of 
Ag. They can be expressed as 


OO 

t{X) = QtX + y''\x) + y^yx) iyy sinriAflAX, 

nA=l 

OO 

(p{X) = n^X + (f^^yx) + if^^yx)] If^^yX) sin ua^aX, 

nA = l 

where the index A runs over {r, 9}, and 



(18) 

(19) 


( 20 ) 


with (•••);, = lim 7 ’_).oo( 2 r)“^ dX - ■ ■, representing the time average along the geodesic. 
We choose the initial conditions as t{X = 0) = if{X = 0) = 0. corresponds to the frequency 
of the orbital rotation. 

In Appendices A and B, we present the PN formulae of the orbital parameters, {E, L, C}, 
the fundamental frequencies, {flj, He; and the Fourier coefficients of the motions in 
Fqs. (15), (16), (18) and (19). 


^ If the ratio of the radial and longitudinal frequencies is irrational, we can adjust the origin of A 
approximately so that the radial and longitudinal oscillations reach the minima simultaneously at 
A = 0 [8]. On the other hand, it is not the case if the ratio is rational, i.e. the resonance case. This 
implies that the secular evolution of a resonant orbit cannot be described only by the PN formulae 
derived in this work [20]. 
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2.2. Secular evolution of orbital parameters 

The gravitational perturbations in Kerr spacetime can be described by the Weyl scalar, T 4 , 
which satisfies the Teukolsky equation [21]. To solve the Teukolsky equation, the method of 
separation of variables is often used, in which '^4 is decomposed in the form as 

= [ duRA{r)SA{e)e^^^-^^\ ( 21 ) 

im 


where Sa{0) is the spin -2 spheroidal harmonics and A represents a set of indices in the 
Fourier-harmonic expansion, {i,m,u}}. The separated equation for the radial function is 
given by 



+ 


/ + 4i(r — 

V ^ 


M)K 


Siujr — A 


RA{r) = Ta, 


( 22 ) 


where 


K = {r"^ + 0 ^) 0 ) — ma, 


Ta is the source term constructed from the energy-momentum tensor of the point particle, 
and A is the eigenvalue determined by the equation for Sa (To find the basic formulae for 
the Teukolsky formalism used in this paper, refer to the section 2 in [22] for example). 

The amplitudes of the partial waves at the horizon and at inhnity are defined by the 
asymptotic forms of the solution of the radial equation as 

= iiA(r ^ 00 ) =(23) 

with r+ = M -|- Vand k = u — mo/(2Mr_|_). Since the spectrum with respect to a; 
gets discrete in the case of a bound orbit, Z^’°° take the form 

= 2n5iuj - (24) 

where A denotes the set of indices, {i,m,nr,ng}, and 

a^mUrne — “f“ T . (25) 


With these amplitudes, the secular changes of the orbital parameters, {E,L,C}, can be 
expressed by 



VE 


VE 


47rcj^ 

^mnr-ne 


m 


^mrirTie 


r/OO 


ryOO 


“1“ n 


+ Oi£jn n 


zf 


A 

) 

2 \ 

^A 



d^\ 

dt /, 


-2{a^Ecos^e)^{^) +2{Lcot^e), 


dt/t 




A 


ne^e 

^TrUJ^n^ng 


ryOO 


T Oilrnip^mnrne') 


z 


H 


(26) 

(27) 


(28) 


where 


Q^£m(^) 


256{2Mr+fk{k‘^ + 4e^){k'^ + 16e‘^)uj^ 


= a/m 2 -aV(4Mr+), 


and Cs is the Starobinsky constant given by [23] 

\^S?‘ = [(A-I- 2 )^-I-datum — [A^-I-36aa;m — 36a^a;^] 


( 29 ) 
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(30) 


+(2A + 3)(96a^a;^ — ASaojm) + 144w^(M^ — a^). 

It should be noted that, in these formulae, the averaged rates of change are expressed with 
respect to the Boyer-Lindquist time, which can be related to those with respect to A [24] as 



for a function of time, I{t). Also it should be noted that each formula in Eqs. (26)-(28) 
can be divided into the infinity part and the horizon part; the former consists of the terms 
including the amplitudes of the partial waves at the infinity, the latter consists of the 
terms including the amplitudes at the horizon, Z?. As for the energy and azimuthal angular 
momentum, the infinity parts are balanced with the corresponding fluxes radiated to infinity 
and the horizon parts with the absorption of the gravitational waves into the central black 
hole [23, 25], 

The practical calculation of Z?’°° involves solving the geodesic equations, calculating two 
independent homogeneous solutions of Eq.(22) and the spin-2 spheroidal harmonics, and the 
Fourier transformation of functions consisting of them. In this work, we followed the same 
procedure proposed in [15] to perform these calculations analytically. 

In performing the summation in Eqs. (26)-(28) practically, we need to truncate the sum¬ 
mation to finite ranges of A = {£,m,nr,ng}. To obtain the accuracy of the 4PN and 0(e®), 
it is necessary to sum i in the range 2 < £ < 6 (2 < ^ < 3), in the range —3 < < 3 

(—2 < Ur < 3) and n$ in the range —8 < ng < 12 (—4 < ng < 6) for the infinity (horizon) 
part. The other modes out of these ranges are the higher PN corrections than the 4PN order 
or the higher order corrections than 0(e®). 


3. Results 

3.1. PN formulae of the secular changes of orbital parameters 

In this work, we derived the analytic 4PN order formulae of Eqs. (26)-(28) in the expansion 
with respect to the orbital eccentricity, e, up to 0(e®) (we simply call them as the 4PN 0(e®) 
formulae). Since the full expressions of the 4PN 0(e®) formulae are too lengthy to show in 
the text, we show the infinity parts up to the 3PN order and the horizon parts up to the 
3.5PN order (while we keep the expansions with respect to e up to 0(e®)). The complete 
expressions of the 4PN 0(e®) formulae will be publicly available online [26]. 

The infinity parts of Eqs.(26)-(28) are given by 


dt ] ^ 


dE 

dt 
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1 + WT e" + — + <^ - 


24 
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+ ^ 47r- —yg + 


96 
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336 


672 
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e^ + 


8609 
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vr- Yq ] e 

48 24 ^ ' 
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3935 949 ^^ A 4 , /10007 491 ^ ^ g I 3 

-TT- Yq\e^-\-{ -vr-To e° 

192 32 V 9216 192 ^ ' 


44711 527, . 2 329 2 f 172157 4379 2 6533 , . 2 \ 2 

“ 9072 ^ { 2592 1^ ^ ^ ® 
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2764345 3823 


24192 


256 


q + 


6753 


Y\^ \ 


7/32 

















+ 


/3743 363 




-^(1 + 


2855 


2304 512 


r 3749 8191 

+ S - 7 :;;;^ VT + ( - 


1536 
44531 




336 


672 


336 


■ VT + 


1759 


Yq\e 


/4311389 111203 

( TT + 


43008 


1344 




56 

/15670391 


49685 


387072 


■ vr — 


448 


Yq e® }v 


+ < 


r 6643739519 1712 


69854400 


105 


7 


3424 , , , 16 p 135 9 169 

- In 2 H-vr^ +- q^ - 

105 ^ ^ 3 8 6 


irYq 


73, . 2 /43072561991 680 o 234009 , , , 14552 

+ — Y^q^ + -H-vr^-In (3)- 

21 ^ V 27941760 9 560 ^ ^ 63 


7 


+ 


13696 , , , 205747 o 4339 13697 ^ . 2\ 2 

-n 2 +-g- 7rYq+ - Y^q^ 

315 ^ ^ 1344 ^ 16 ^192 ^ 

, 2106081 , , , 12295049 , , , 

vr" +- n (3-n (2 

V ; ^260 ^ ’ 


/919773569303 5171 


\ 279417600 


+ 


36 


553297 


■7 + 


208571 


448 
42271 


q 


1260 ' ' 1792 " 96 

/308822406727 864819261 


471711 . 2 

TTYq-\ --- Y^q^ \e‘ 


1792 


+ 


V 186278400 

5224609375 

193536 

289063 ,,-^ 

H- Y^q 

1536 


35840 


In (3) - 


187357 


, , , 24908851 , 

In (5) +-—-- In (2) + 


1260 
3253 


7 + 


1751 


36 


TT 


252 


10752 


4867907 

27648 


■irYq 


2.2 


1712 14552 2 553297 . 187357 . > , . 

+ ^ 77 ^ + ..... + __ e® 1 Inn S>n® 


105 


63 


1260 


1260 


di /n 


, 7 nK. f 1247 425 n 10751 4 i o 

l + -e^}Y+< -H- e^\Yv^ 

8 I 1 336 336 2688 




63 97 91 ,,2 \ 2 

— q^ -Try- Y^q 

8 8 4 ^ ' 


95 


49 


461 


49 


vrye® 


+ + V" - 4608 

44711 57 2 45 , . 9 f 302893 201 2 37, . 

+ <! -TTTTz:^ - 7779 + TT ^ q + -^777-T7r9 + 77 Y^Y 


9072 16 


8 


6048 


16 


+ < 


n- 

r 4301 


701675 351 2 331 2 2\ 4 162661 


224 


- q^ + - Y^q^ 

24192 128 64 ^ 

,.2 8191 2633 

y 9 - t^ttY - 


e" + 


16128 


yn^ 


672 


224 


, 66139 

’ - 1344 ^ 


48361 18419,0 \ 2 

■ Try H- tttt- Y^q 


1344 


'3959 1657493 

+ I T777779+ _ Try- 


+ 


1792 


/19161 


I 


3584 


9 + 


43008 
5458969 
774144 


448 

257605 

5376 


Y\ ] 


52099 2 \ 6 I 5 

Try- Y^q e® }v^ 

1536 ' ‘ 


(32) 


8/32 

































































145 


6643739519 


S - TT O H” 

^ 12 ^ 69854400 


■ y + —TT^y- 


1712 

105 


7y - 


3424 


In ( 2 ) y 


171 


145 


1769 


- Yq^ - 7rY^q + 

112 ^ 4 ^ 112 


3^2 




/995 229 2 ,. 6769212511 1391 , , ,,, 24503 

+ ( — vr g + — vr^y + y + In (2) Y - 7 Y 


12 
78003 


6 


8731800 


30 


210 


, , ,,, 46867 , o 877 27997 ,0 2 \ 2 

- n 3 y- Yq^ -tt Y^q +- Y^q^ 

280 ^ ’ 1344 ^ 4 ^192 ^ 


/ 21947 4795392143 3042117 , , , 418049 , , , 

+ 1 -^T 7 T^ 7 rg+ y+ ln(3)y-- In (2) y 


V 384 
11663 


7761600 


1120 


84 


109 2 ,. 1481 2 22403 , . 267563 ,,3 2 \ 4 

- 7 y H-TT^y-yg-vry^gn-y'^g^ 

140 ' 4 16 ^ 128 ^ 1344 ^ 

/38747 31707715321 23 2 .. 94138279 , , ,,, 

+ -TT g +-y H-TT^y +- In 2) Y 

1 13824 ^ 10^0^70.nn ^ o-icn W 


186278400 


16 


2160 


1044921875^ 42667641 , , 2461 68333 ^^2 

in(5)y--— in(3)y- 7 y- ^ 7777 -yg 


96768 


3584 


560 


3584 


59507 , . 183909 , o 2 

■ ^ ^ <1 + oro. ^ Q I 


4608 " 3584 

/1712 24503 o 11663 


■ 


105 


+ 


210 


-e" + 


dt 


N 


1 + 7 + 1 - 

85 


+ <j47r —7-yg + ( —TT- — yg ) e 

49 


140 

743 23 

^ 42 

97 211 


+ 


2461 


560 


e® I yinv U® 


(33) 


+ 776 ^ + 


11927 

2688 


4 \ 2 

e V 


7 6® 


V32 64 4608 

129193 329 2 53 , . 2 f 84035 929 2 163, . 2 \ 2 


+ - 


1030273 1051 


48384 


2 , 387 ^2 2 \ 4 , 100103 6 I 4 

-g^ H- Y^q^ H- )-v^ 

768 ^ 64 ^ / 8064 


+ — 


4159 

^12 


TT + 


2553 


Yq+\- 


21229 


553 


■ 7 - Yq\ e 

1344 192 ^ ' 


/2017013 475541 \ 4 /6039325 153511 \ « 1 ^ 

+ 1 7 —yg 6^ + ( 7- yg 6° >f® 


I- 


43008 


5376 


11683501663 16 2 1712 

^ 139708800 3 105 


V 774144 
3424 


3584 


193 


7- 


In (2) + ^q^-—7rYq 

105 ^ ^ 192 ^ 4 ^ 


+ - 


2515 


y^g^ + I 


48 


/16319179321 ^ 229 _2 24503 
23284800 ^ 210 


1391 , , , 
7 H-In ( 2 ) 

' 30 ^ ^ 


78003 

280 


In (3) + 


/ 211889615389 109 2 

Y I -4“- 7 “l" 

V 372556800 


16979 

1344 

LOG 

T 


2077 118341 , 0 2 \ 2 

7yg H-—— y^g^ 16^ 


8 

3042117 

1120 


448 




140 


84 


9/32 


































































132193 2 24543 91747 , ^ 2 \ 4 

" ■T^Yq+ Y^q^]e^ 


3584 


128 


336 


/ 33928992071 1044921875 , , , 23 2 42667641 , ^ ^ 

+ V 186278400 - 9676^+ 16 ' “^5^'"® 

94138279 , , , 2461 24505 2 4151 718799 ,0 2 \ r 

+ In (2) - 7 - q - ^777 '^Yq+ -77777 Y^q^ ) e® 


2160 ^ ^ 560 

/1712 24503 o 11663 


5376 


288 


10752 


- 


+ 


■e^ + 


105 ' 210 ' 140 

where the leading contributions are given by 

'dE\ 


■e^ + 


2461 


\\d.v\v 


,,6 


(34) 




dt J 


f fA)%10(l_e2)3/2^ 


N 


5 \MJ 


dC 


5 \M 


dt 


64 




(35) 


N 


The horizon parts of Eqs.(26)-(28) are given by 
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{dE/dt)f^, {dL/dt)Y and {dC/dt)^ in Eqs. (36)-(38) are new PN formulae derived in this 
paper. (dE/dt)'^, {dL/dt)^ and (dC/dt)^ in Eqs. (32)-(34) are consistent with those in 
Ref. [15] up to 2.5PN and O(e^). 

From the leading order expressions in Eq. (35), one will find that the Carter parameter, 
C, does not change due to the radiation of the gravitational waves when Y = 1 (equatorial 
orbits) because (dC/dt)^ = 0. 

In the Schwarzschild case, the Carter parameter corresponds to the square of the equato¬ 
rial angular momentum (the normal component to the rotational axis of the central black 
hole). Then there is expected to exist the duality between and C due to the spherical 
symmetry. In fact, from Eqs. (33) and (34), (and also from (37) and (38)), one can find that 
{dL?'/dt)t vanishes in E = 0 (polar orbits) while {dC/dt)t for Y = 0 coincides with {dL‘^/dt)t 
for y = 1. This can be also realized by seeing that the secular change of the total angular 
momentum, {d{L?‘ + C)/dt)t, is independent of Y . Then, it might be possible to understand 
that {dL/dt)f‘ becomes 1.5PN from the leading order when ^ 7 ^ 0 and E = 0 (polar orbits) 
due to the spin-orbit coupling. 

From the expressions of the horizon parts shown in Eqs. (36)-(38), we find that the absorp¬ 
tion of the gravitational waves to the central black hole contributes at O(v^) from the 
leading order in Eq. (35) for g / 0 and at O(v^) for q = 0. The O(v^) and 0(v^) corrections 
in {dE/dt)Y can be positive for <? > 0 , which means that the particle can gain the energy 
through a superradiance phenomenon. These observations are consistent with the results for 
circular, equatorial orbits shown in Refs. [27, 28]. 

We also find that the superradiance terms in Eq. (36) vanish for E = 0, and that {dE/dt)^ 
has only the 4PN and higher order corrections. The superradiance terms may come from 
the coupling between the black hole spin and the orbital angular momentum, like oc L • S oc 
qcosL. Hence, when the orbital inclination increases (E gets small correspondingly), the 
superradiance is suppressed [29]. 

3.2. Comparison to numerical results 

To investigate the accuracy of the 4PN 0(e®) formulae derived in this work, we compare 
them to the corresponding numerical results given by the method established in Ref. [16-18], 
which enables one to compute the modal fluxes with the relative error of 10 in double 
precision computations. In the practical computations, as well as in deriving the analytic 
expressions, we need to truncate the summation to finite ranges of A = {i,m,nr,ng} in 
Eqs. (26)-(28). In order to save the computation time in the numerical calculation, we sum 
I up to 7. We can check that the error due to neglecting terms for £ > 8 is smaller than 
the relative error in the 4PN 0(e®) formulae from the corresponding numerical results up 
to £ = 7. We also truncate Ur and ng to achieve the relative error of ~ 10“^ in numerical 
results up to £ = 7. For the parameters investigated in the comparison, the relative error 
of ~ 10 “'^ achieved by truncating and ng is again smaller than the relative error in the 
4PN 0(e®) formulae from the numerical results up to ^ = 7. Thus, we can regard numerical 
results as benchmarks to investigate the accuracy in our analytic formulae. 


11/32 




Here we define the relative error in the analytic formula of {dE/dt)t by 


A.E = 


1 - 


\ dt A / \dt A 


(39) 


where {dE/dt)^^^ denotes the analytic formula in order to distinguish it from the corre¬ 
sponding numerical result, We also define the relative errors in the analytic 

formulae of {dL/dt)t and {dC/dt)t in a similar manner and denote them as and Ac 
respectively. 

Fig. 1 shows several plots of A^ for the 4PN 0(e®) formula as a function of p for several 
sets of (e, i) with q = 0.9. In the plots, we also show the relative errors in the 2.5PN O(e^) 
and 3PN O(e^) formulae for reference. From the plots for e = 0.1 (three on the top), one can 
find that A^ falls off faster than for p > 10 (Similarly, the relative errors in the 2.5PN 
O(e^) and 3PN 0{e^) formulae fall off faster than and p~^). Noting v = \/T/p, this 

would be a good indication that our PN formula has been derived correctly up to required 
order. 

Ac is expected to contain not only higher order corrections than the 4PN order, but also 
the higher order corrections of eccentricity than 0(e®) in the lower PN terms, which will 
become dominant when p and e get larger. In fact, seeing the plots for e = 0.7 in Fig. 1, 
one can find that the relative error strays out of the expected power law line for large p. 
This behavior is clearer in the plots of the relative error in the 2.5PN O(e^) formula. From 
Eqs. (32) and (36), we know that the relative error in the 2.5PN O(e^) formula contains 
the 0{e^) correction in the 0(u®) term. The effect of this correction appears as large-p 
plateaus in the plots (also see Fig. 6). This may motivate us to perform the higher order 
expansion with respect to the orbital eccentricity in the PN formulae or to derive the PN 
formulae without performing the expansion with respect to the orbital eccentricity [30-33]. 
In addition, it might be noted that the behavior of the relative error does not strongly 
depend on the inclination angle i for fixed q and e. 

In Fig. 2, we show the relative errors in the 4PN O(e^) formulae for the secular changes of 
the three orbital parameters, {E,L,C}, for several sets of {q,e) and l = 50°. As in the case 
of Ag shown in Fig. 1, the relative errors. Ax, and Ac, fall off faster than p~‘^ when p > 10, 
except for the large p region [p > 100) in the case of e = 0.7. Thus, the 4PN 0(e®) formulae 
for the secular changes of the orbital parameters are expected to be valid up to 0(u®). From 
Fig. 2, one might think that it is enough to investigate only Ag to discuss the accuracy of 
our formulae since there are not large differences in the relative errors, Ag, Ax and Ac- 

Fig. 3 shows contour plots for Ax as a function of p and e for several sets of (t, q). From 
these plots, one may be able to comprehend the accuracy of our PN formulae more easily 
than using Figs. 1 and 2. One will find that the relative error becomes smaller (larger) for 
larger (smaller) p and smaller (larger) e. Moreover, it might be noticed that the relative 
error does not strongly depend on the inclination angle l for fixed q as expected from Fig. 1. 
If one requires Ax < 10“^ as an error tolerance, one can use the contour line with the label 
10“^ to find the region of validity in the figure. For example, one will find that Ax < 10“^ 
for p > 50 and e = 0.1, p > 80 and e = 0.4, and p > 120 and e = 0.7 when q = 0.9. 
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Fig. 1 The relative errors in the analytic PN formulae for the secular change of the 
particle’s energy as a function of the semi-latus rectum p for g = 0.9, e = 0.1,0.4 and 0.7 
(from top to bottom) and i = 20°, 50° and 80° (from left to right). In addition to the error 
in the 4PN 0(e®) formula, those in the 2.5PN O(e^) and the 3PN O(e^) formulae are shown 
in each plot for reference. We truncated the plots at p = 6 because the relative errors get 
too large (nearly or more than unity) in p < 6 to be meaningful. One finds that the relative 
error becomes smaller with increasing orders of the PN approximation and the expansion 
with respect to the eccentricity. The relative error in our 4PN 0(e®) formula falls off faster 
than p“^ when the eccentricity is small, e.g. e < 0.4. Since v = this would imply that 

our 4PN formula is correctly representing the secular change up to the 4PN order. Note, 
however, that the relative error in the 4PN 0(e®) formula for e = 0.7 falls off slower than 
p~^ when the semi-latus rectum becomes larger,e.p. p > 100. This might be because of the 
higher order corrections of e than 0(e®), which will contain the lower PN terms than the 4PN 
order. We also note that changing the inclination angle, i, does not change the dependence 
on p of the relative error for hxed q and e so much. This might be checked more easily in 
contour plots in Fig. 3, which show the relative error as a function of p and e for fixed q and 

i. 


3.3. Implementation of an exponential resummation method 

In order to improve the accuracy in the analytic PN formulae, one may apply some resum¬ 
mation methods such as Pade approximation [34], the factorized resummation [35-37] and 
the exponential resummation [38]. Since the exponential resummation may be the simplest 
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Fig. 2 The relative errors in the analytic PN formulae for the secular changes of the 
three orbital parameters, {E, L, C}, as functions of the semi-latus rectum p for l = 50°, q = 
0.9, 0.5,0.1 and —0.9 (from top to bottom) and e = 0.1,0.4 and 0.7 (from left to right). We 
truncated the plots at p = max{6,ps(e, i)}, where Ps(e, t) is the value of p at the “separatrix” 
(the boundary between stable and unstable orbits), because the relative errors get too large 
in p < 6 to be meaningful and the orbit is not stable for p < Ps(e, t). As pointed out in Fig. 1, 
the relative errors in our 4PN 0(e®) formulae fall off faster than p~^ when the eccentricity 
is small, e.g. e < 0.4, while the fall-off gets slower when p is larger for e = 0.7. There are not 
large differences in the behaviors of Ag, and A^. This suggests that it might be enough 
to focus only on {dE/dt)t to investigate the accuracy and convergence of our 4PN formulae. 


one to implement among them, we here choose to implement the exponential resummation. 
We apply it to our 4PN formulae and check how the accuracy is improved. 

To introduce the exponential resummation, we make use of the following identity 
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Ae lor q=0.9 and i=20" 


Ag lor q=0.9 and i=50" 


Ag for q=0.9 and 1=80° 



Fig. 3 The relative error in the 4PN 0(e®) formula for the secular change of the particle’s 
energy, A^;, as a function of the semi-latus rectump and the eccentricity e for q = 0.9,0.5, 0.1 
and —0.9 (from top to bottom) and i = 20°, 50° and 80° (from left to right). We truncated 
the plots at p = max{6,ps(e, t)} because the relative errors get too large in p < 6 to be 
meaningful and the orbit is not stable for p < ps[e,i). From the figures, it is easily found 
that the relative error becomes smaller (larger) for larger (smaller) p and smaller (larger) e 
with fixed q and t. If one requires the relative error to be less than 10“^, the region in the 
semi-latus rectum p and the eccentricity e will be p > 50 and e = 0.1, p ^ 80 and e = 0.4, 
and p > 120 and e = 0.7 when q = 0.9. It might be noticed that the relative error does not 
strongly depend on the inclination angle i for fixed q as pointed out in Fig. 1. 


where / = {F,L,C'}. The exponential resummation can be obtained by replacing the 
exponent in (40) to the expansion with respect to u, 




truncated after nth order of v 


(41) 
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where we do not perform the expansion with respect to e. Since our PN formulae for {dl/dt)t 
are given at the 4PN order, we truncate after 0(n®). Finally, the exponential resnmmed 
form is expressed as 


dt/t 


dt 


exp F/ 


(42) 


N 


Fig. 4 shows the relative errors in the exponential resummed forms of the secular changes 
of E, L and C, estimated by using Eq. (39). We also show the relative errors in the Taylor- 
type formulae in the same graphs for comparison. One will find that the relative errors in 
the exponential resummed forms are less than those in the Taylor-type formulae in most 
cases, except for {dC/dt)t in the case of g = 0.9, (e, t) = (0.1,50°). Using the exponential 
resummation when q = 0.9 and i = 50°, the region to satisfy A^; < 10“® is extended to 
p > 40 from p > 50 for e = 0.1, p > 60 from p > 80 for e = 0.4, and p > 100 from p > 120 
for e = 0.7. This might motivate ns to use the resummation method to improve the accuracy 
of Taylor-type formulae even in the case of general orbits. 


3.4- Convergence with respect to v and e of the analytic formulae 
Apart from comparisons to the nnmerical results, we may also discuss the convergence prop¬ 
erty in our PN formulae with respect to v and e by investigating the contribution of each 
order of v and e in the formulae although this is a rough estimation. 

First we assess the PN convergence of onr formulae. For this purpose, we introduce as 



where p=l/n^ and A„ is the 0(u”') term in the PN formula of {dE/dt)t, e.g. Aq = 
1 + H -I- II Ai = 0 and A 2 = + fit + UtI depends on 

{q,p, e, Y) in general although we omit the argument for simplicity. 

Since A„ shows the relative importance of the 0(n”) term in the PN formnlae, it can be 
nsed to investigate the convergence property with respect to v: it is expected that |A„_|_i| < 
|A„| for moderately large n if the PN formnla converges. In Fig. 5, we plot the relative 
contribution of each order, A„, as a fnnction of p for several sets of (e,t) and q = 0.9. 
From this fignre, one may find that A„ does not strongly depend on the inclination angle, 
i, as shown in Sec. 3.2, while it strongly depends on e. The convergence gets worse when 
the orbital eccentricity becomes larger. This tendency is particularly evident in the small-p 
region. Fixing the value of p, the orbit with larger e passes by closer to the central black 
hole and will be affected by the stronger gravitational field. Hence the PN convergence is 
expected to be worse when the eccentricity becomes larger. 

Next, in order to investigate the convergence of the expansion with respect to the orbital 
eccentricity in the PN formula, we introduce An as 

where the term Aq coincides with the energy flnx for circular orbits and An = 0 when n is 
odd. 

One may ask whether the condition, \A 2 n+ 2 e‘^'^~^‘^\ < |^ 2 ne^”|, is satisfied for moderately 
large integer n if the series converges. From Fig. 6, it is found that the condition is satisfied 
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Fig. 4 The relative errors in the PN formulae and the exponential resummation formulae 
for the secular changes of the orbital parameters, {E,L,C} as functions of the semi-latus 
rectum p ioi q = 0.9, l = 50° and e = 0.1, 0.4 and 0.7 (from top to bottom). We truncated 
the plots at p = 6 because the relative errors in the PN formulae get too large in p < 6 to 
be meaningful. Using the exponential resummation, the accuracy is improved in most cases. 
For example, the region to satisfy Ag < 10“^ is improved from p ^ 50 to p > 40 for e = 0.1, 
p > 80 to p > 60 for e = 0.4, and p > 120 to p > 100 for e = 0.7. This would suggest us to 
try to apply resummation methods to the PN formulae even in the case of general orbits. 

in most cases. As expected, the convergence becomes slower when the eccentricity is larger. 
Especially, the convergence gets worse when p < 10 in e = 0.7 case. The calculation of the 
higher PN corrections will be necessary to improve the bad convergence for small p. We also 
note that An does not strongly depend on l for a fixed q as in Sec. 3.2. 

4. Summary 

We have derived the secular changes of the orbital parameters, the energy, azimuthal angular 
momentum, and Carter parameter of a point particle orbiting a Kerr black hole, by using the 
post-Newtonian approximation in the first order black hole perturbation theory. We have 
extended the previous work [15], which derived formulae up to the 2.5PN order with the 
second order correction with respect to the eccentricity, to the 4PN order with the sixth 
order correction with respect to the eccentricity. We have also included the contribution 
due to the black hole absorption, which has not been included in [15]. As shown in the 
case of equatorial, circular orbits [27, 28], we have found that the secular changes of the 
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Fig. 5 The relative contribution of the 0(n") term in the PN formula for {dE/dt)t, defined 
in Eq. (43). We plot the absolute value of as a function of the semi-latus rectum p for 
e = 0.1,0.4 and 0.7 (from left to right), and i = 20°,50° and 80° (from top to bottom) 
when q = 0.9. We truncated the plots at p = 6 because the relative contributions get too 
large in p < 6 to be meaningful. It is expected that |A„+i| < |A„| for moderately large n 
if the PN formula converges. As shown in Sec. 3.2, A„ does not strongly depend on i for a 
fixed e although A„ strongly depends on e. In fact, the convergence seems worse the orbital 
eccentricity becomes larger. This tendency is clear for small p, e.g. p < 10. 


three orbital parameters due to the absorption, appear at the 2.5PN (4PN) from the leading 
order in the Kerr (Schwarzschild) case, and that the 2.5PN and 3.5PN contributions of the 
absorption to the secular change of the particle’s energy can be positive for q > 0, which 
implies that a superradiance can be realized in the Kerr case. We have also found that 
the superradiant contributions in the secular change of the energy get smaller when the 
inclination angle becomes larger and they vanishes for polar (Y = 0) orbits. This means 
that the superradiant scattering may be suppressed for inclined orbits [29]. 

To investigate the accuracy in our 4PN formulae, we have compared the formulae to high- 
precision numerical results [18] in Sec. 3.2. We have found that the accuracy gets worse when 
the orbital velocity and the orbital eccentricity become larger, as expected. If the relative 
error in the 4PN 0(e®) formula for the secular change of the energy is required to be less 
than 10“^, the parameter region to satisfy it might be p > 50 for e = 0.1, p ^ 80 for e = 0.4, 
and p > 120 for e = 0.7 when q = 0.9. The region does not strongly depend on the orbital 
inclination angle. From Fig. 1, one can clearly find the improvement of the accuracy in our 
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Fig. 6 The relative contribution of the 0(e"') term in the PN formula for {dE/dt)t, defined 
in Eq. (44). We plot the absolute value of AnC^ as a function of the semi-latus rectum p 
for L = 20°, 50° and 80° (from left to right) when q = 0.9. We truncated the plots at p = 6 
because the relative contributions get too large in p < 6 to be meaningful. It is expected that 
|^ 2 n+ 2 e^"'"’~^| < |^ 2 ne^”| for moderately large n if the series with respect to e converges. This 
condition is satisfied in most cases shown in this figure. The convergence becomes slower 
when the eccentricity is larger. Especially, the convergence for p < 10 is quite bad in e = 0.7 
case. We also note that An does not strongly depend on i for a fixed q as mentioned in 
Sec. 3.2. 


PN formulae from the previous work at the 2.5PN order and the second order correction in 
the orbital eccentricity [15] whose relative error is larger than 10“^ for p > 100 when e > 0.4 
since, in this region, the error due to the truncation of the expansion with respect to the 
orbital eccentricity is larger than the one of the PN expansion. 

One may improve the accuracy of our PN formulae by using resummation methods. In 
this paper, we have applied the exponential resummation [38] to our 4PN formulae and 
confirmed that the resummation method improves the accuracy in most cases investigated 
here. For example, we found that the region in which the relative errors are less than 10“^ 
can be extended from p > 50 to p > 40 for e = 0.1, p > 80 to p > 60 for e = 0.4, and p > 120 
to p > 100 for e = 0.7. 

We also investigate the convergence properties of the PN expansion and the expansion 
with respect to the orbital eccentricity, respectively. Both convergences get worse when the 
semi-latus rectum is smaller; in other words, the gravitational field becomes stronger. This 
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tendency gets clearer in the case of large eccentricity, in which the particle passes by closer 
to the central black hole. 

In order to improve the accuracy and convergence of the 4PN 0(e®) formulae near the 
central black hole and to obtain the physical information of the source in the strong-held 
region, it is necessary to derive the higher order corrections of the PN expansion and the 
expansion with respect to the eccentricity. It may be possible to avoid the expansion with 
respect to the eccentricity and to derive the PN formulae applicable to arbitrary eccentricity. 
So far the PN formulae of the rate of the energy loss without performing the expansion with 
respect to the eccentricity had been derived for equatorial orbits in [30-33]. The extension of 
these results to the case of inclined orbits is challenging: we can obtain analytic expressions 
for general bound geodesic orbits in Kerr spacetime without performing the expansion with 
respect to the eccentricity nor the inclination by using results in Ref. [39], while we need 
to reformulate the source term of the Teukolsky equation and the derivation of the partial 
waves constructed form the source term. We would like to leave it to the future work. 
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A. PN formulae for the orbital parameters and fundamental frequencies 

In this section, we present the PN formulae of the orbital parameters, {E,L,C}, and the 
fundamental frequencies, Here we show the formulae up to the 3PN 0(e®) 

order to save space although it is possible to calculate them to the higher order. The higher 
order results will be publicly available online [26]. 
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B. Fourier coefficients of bound orbits 

Here we show the PN formulae of the Fourier coefficients in Eqs. (15), (16), (18) and (19) 
up to the 3PN 0(e®) order. The 4PN 0(e®) results obtained in this work will be available 
online [26]. 

In this work, we follow the same procedure as in [15] to derive the amplitudes of the partial 
waves, Z^’°° in (23). In the formal expression, the dependence of appears in the form 
of the combination as X = sin ^e**^^**', which can be expressed in the Fourier series as 


00 

X = p'^ cos ngQgX + iX^g sin ngXlg^ , 

ns=0 


(Bl) 


Therefore we show the Fourier coefficients of X instead of 
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B.2. Longitudinal component 
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B.3. r-part of the temporal component 
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B.Ji^. 9-part of the temporal component 


— t 
P 


m _ 


= 0 , 


1 fW _ 
''2 ~ 
P 

1 


fW = 


P 


(y2 -1) 2 3 (y2 -1) 2 5 

1 --- -q^V^ - - - -q-^v^ 


0 (n: odd) 

(^(^2n-l) gygj^^ 


y(y2 - 1) 


(3 + e'^)q^v^ 


(B22) 

(B23) 

(B24) 


B.5. r-part of the azimuthal component 
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B.6. 9-part of the azimuthal component 
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C. Secular evolution of the orbital parameters, v, e, and Y 

An alternative set of the orbital parameters, J = {v, e, Y }, is also useful to specify the orbit. 
The secular changes of the parameters can be derived from those of / = {E, L, C}, as 



E (G-‘)/ 

I=E,L,C 



(Cl) 
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where Gj = d{E,L,C)/d{v,e,Y) is the Jacobian matrix for the transformation from 
{E,L,C} to {v,e,Y} 

Substituting the 3PN 0(e®) formulae of {dl/dt)f^ shown in Sec. 3 into the above relation, 
we obtain the secular changes of {v, e,Y} associated with the flux of gravitational waves to 
inhnity as 
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^It should be noted that, to calculate the Jacobian matrix up to 0(e®), one need to calculate 
{E,L,C} up to 0(e®) since the leading terms do not depend on e and then the relative orders of 
accuracy of their first derivatives with the eccentricity of the Jacobian matrix in Eq. (Cl) are reduced 
by O(e^). For a similar reason, one also need to calculate E up to 5PN order because the relative 
PN order of dE/dv is reduced by 0{v^) compared to E. 
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where the leading contributions are given by 
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In the same way, substituting the 3.5PN 0(e®) formulae of {dI/dt)Y shown in Sec. 3 into 
Eq. (Cl), we obtain the secular changes of {u, e,Y} associated with the flux of gravitational 
waves to the horizon as 
H 


dt/t 


dv 

dt 


N L 


-^{8 + 24e2 + 3e^}{8 + 9g2 + 15y2g2| 

256 


11 189 2 15 ^2 2 


64 


64 


69 81 2 45^0 


9- 777 9 + -^-7r 9+7777 Y^q 


381 7479 2 1035 , ^ 2 


, 11 423 2 45 ^ . 2 

+-rn- +9 

' 32 512 ^ 512 ^ 


32 


Y qv' 


(C6) 


de \ 

dt/t 


H 


de\ 

dt). 


33 

4864 


{8 + 12 e 2 + e ^}{8 + 952 y i 5 y 2 ^ 2 | y^^s 


453 8127 ^2 45 

152 ~ 1216 ^ ~ 1216 


Y^q^ + ( - 


2979 28593 


304 


1216 


9 + 


1485 


2„2 


Y^q 


5649 111591 2 16515 2 2 

+ 1 -T77777 - ^777777^ 9 + 77777^ + 9 


dt 


H 




1216 9728 

3 


9728 


Yqv' 


(C7) 


N L 


7808 


{8 + 24+ 3} {16 + 33+ 15Y^q^j 


51 585 , o 2 1953 

-_ 1 _-y n —- 

122 1952 ^ 1952 


225 16875 2 3375 2 2 


^ ' 61 1952 


9 " + 


1952 


Y^q^ 


+ - 


2961 109863 2 16335 . 2A 4 


976 


15616 


9 + 


15616 


Y^q^ 


171 3159 2 405 , . 2 

+ 1 -777777 - 77777777 9 + 77777777 + 9 


976 7808 


7808 


(C8) 


Actually, we can obtain the higher PN results by using the 4PN O(e^) formulae for the secular 
changes of {J^,L,C}, although we do not present them in the text. The full expressions of 
{dJ/dt)^ and {dJ/dt)^ for J = {v,e,Y} will be available online [26]. 

Here we make a comment on the reliable order of the expansion with respect to e in 
{de/dt)t- By using Eq. (Cl), {de/dt)t can be calculated from the linear combination of 
the secular changes of {E,L,C}. Since the leading order of the (e,/)-component of the 
inverse Jacobian matrix is 0(l/e), each term in the linear combination is apparently 0(l/e). 
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Fig. Cl The relative errors in the analytic formulae for the secular change of the orbital 
eccentricity due to the gravitational waves to infinity. We plot Ag, defined in a similar manner 
to Eq. (39), as a function of the semi-latus rectum p for q = 0.9, e = 0.1,0.4 and 0.7 (from 
left to right) and i = 50°. We truncated the plots at p = 6 because the relative errors get 
too large in p < 6 to be meaningful. The relative error in the previous 2.5PN O(e^) formula 
given in [15] strays off the p~^ line earlier than the 2.5PN O(e^) formula in this paper. This 
trend is clearer for larger e. The relative errors in the 3PN O(e^) and 4PN O(e^) formulae 
fall off faster than p~^ and for small e cases as expected, while this is not the case for 
e = 0.7 because of the higher order correction of e than 0{e^). 


However, the 0(l/e) contribution turns out to vanish due to a cancellation in taking the 
combination, and hence {de/dt)t is 0(e), which corresponds to the well-known fact that 
circular orbits remain circular [40, 41], i.e. {de/dt)t = 0 when e = 0. This cancellation reduces 
the reliable order in {de/dt)t by O(e^), compared to the order of {dl/dt)t for I = {E, L, C}. 
Since we calculate {dl/dt)t up to 0(e®) in this paper, we can obtain {de/dt)t correctly up 
to 0{e'^) from the leading order. 

(dv/dt)^ and (dY/dt)^ in Eqs. (C2) and (C4) are consistent up to the 2.5PN O(e^) order 
with the previous results in Ref. [15], while we find inconsistency in the O(e^) terms of the 
formula for {de/dt)f^ in [15]. This may be explained by the reduction in the reliable order 
mentioned above: the calculations of {dl/dt)t in [15] are done up to O(e^), and therefore 
the resultant formula of {de/dt)t is reliable only at the leading order. We can also confirm it 
numerically. In Eig. Cl we show the relative errors in the two analytic formulae by comparing 
to numerical results [18] in a similar manner to Eq. (39). It can be found that the relative 
error in the previous 2.5PN O(e^) formula strays out of the expected power law, earlier 
than that in our 2.5PN O(e^) formula. This trend is clearer for larger eccentricity. We can 
also confirm the validity of our formula by seeing that the relative errors in our 4PN 0{e^) 
formula falls off faster than p~^. (The leading PN order of the difference in the O(e^) terms 
between the previous and our formulae is If our formula contains any error in the 

term, the relative error will not fall off faster than p~^ for large p.) 

A similar reduction in the PN order occurs in the calculation of {dY/dt)f^: although each 
term in the linear combination of Eq. (Cl) for J = Y is 0{v^), the terms at the first two 
orders, 0(u®) and 0(u^®), vanish due to a cancellation in taking the combination. As a 
result, the leading order of {dY/dt)t is 0(u^^) and hence the reliable order relative to the 
leading term is reduced to 0{v^) (2.5PN order) when we have {dl/dt)t for I = {E,L,C} up 
to 0(u®) (4PN order). 
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Fig. C2 The relative errors in the 4PN formulae for the secular changes of the orbital 
parameters, {v, e, Y }. We plot the relative errors, A j defined in a similar manner to Eq. (39), 
as functions of the semi-latus rectum p for q = 0.9, e = 0.1,0.4 and 0.7 (from left to right) 
and i = 50°. We truncated the plots at p = 6 because the relative errors get too large in 
p < 6 to be meaningful. and Ag fall off faster than while Ay approximately fall off 
as p~^, slower than 0{p~^). This confirms that the relative order of the PN correction of 
the analytic formula for {dY/dt)t is reduced from 4PN to 2.5PN because of the cancellation 
of the low PN terms. 


In Fig. C2, we show the relative errors in the analytic PN formulae for the secular changes 
of the orbital parameters, {v,e,Y}, derived from the 4PN 0(e®) formulae of {dl/dt)t for 
I = {E, L, C}. Similarly in Fig. 1, the relative errors in the analytic formulae for {dv/dt)t and 
{de/dt)t as functions of the semi-latus rectum p fall off faster than p~‘^ when the eccentricity 
is small. Observe, however, that the relative error in the analytic formula for {dY/dt)t falls 
off faster than p~^^‘^, but slower than p~^. 

From the leading order expressions in Eq. (C5), one will find the well known fact that 
equatorial orbits stay in the equatorial plane [10, 15, 42], i.e. {dY/dt)t = 0 when Y = 1. 
In the Schwarzschild case (g = 0), the secular changes of v and e do not depend on Y in 
addition to {dY/dt)t = 0. This implies that the orbital plane can be fixed on the equatorial 
plane {9 = vr/2) due to the spherical symmetry of Schwarzschild spacetime. 

One will also find that the radiation reaction reduces the orbital eccentricity and increases 
the orbital velocity since {de/dt)-^ < 0 and {dv/dt)^ > 0 [40, 41], while the radiation reaction 
increases (decreases) the inclination angle since (dY/dt)^ < 0 {{dY/dt)^ > 0) when q >0 
{q < 0) [10, 15, 42]. Moreover, the secular change of the inclination angle is smaller than those 
of the other orbital parameters since (dlne/dlnu)t = 0{v^) and {dlnY/dliiiv)t = O(v^). 
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